
function [T,x,P,pi,Lp,Ls,gammaL,gammarn,mun,gammaL1,gamma4,sextax, fval] = eval_ms_mc_gov_ss2_nogL(xval)



% for multi-sectors
global rbest  L theta rbest ns n alpha betas d sigma beta ilog_rL delta

% guess wage and share of labor in each sector

w(1:n,1) = xval(1:n); % row vector


for i=1:n
    if (ilog_rL==1)
        rL(i,1:ns) = L(i)*exp(xval(n+(i-1)*ns+1:n+i*ns))./(1+exp(xval(n+(i-1)*ns+1:n+i*ns)));
    else
        rL(i,1:ns) = L(i)*xval(n+(i-1)*ns+1:n+i*ns);
    end
end


T = zeros(n,ns);
for i=1:n
    T(i,1:ns) = alpha(i,:).*rL(i,:)/delta;
end

sextax = zeros(n,n,ns);
for i=1:n-1
    sextax(i+1,1,:) = xval(n+n*ns+(i-1)*ns+1:n+n*ns+i*ns);
end


staf = zeros(n,n,ns);
for i=1:n-1
    staf(1,i+1,:) = xval(n+(2*n-1)*ns+(i-1)*ns+1:n+(2*n-1)*ns+i*ns);
end

pi = zeros(n,n,ns); % (1,2,3): country 2 export to ecountry 1 in sector 3, NOTE country order is opposite to our note
%d, matrix of n times n, dij country j ship to country i
for iss = 1:ns
    
    xnimatrix1= T(:,iss)'.*((w'.*(1+squeeze(sextax(:,:,iss))).*(1+squeeze(staf(:,:,iss))).*d).^(-theta));
    xn1 =sum(xnimatrix1,2).*ones(n,n);
    
    pimatrix1 = xnimatrix1./xn1;
    pi(:,:,iss) = pimatrix1;
    xnimatrix(:,:,iss) = xnimatrix1;
end

betas_ns = ones(n,1)*betas;

x  =  (1+theta)/theta*w.*(L-sum(rL,2));
extax1 = squeeze(sextax(:,1,:));
for m=2:n
    for j=1:ns
        x(1)=x(1)+betas(j)*extax1(m,j)/(1+extax1(m,j))*pi(m,1,j)*x(m);
    end
end
tempt=0;
taf1 = squeeze(staf(1,:,:));
for m=2:n
    for j=1:ns
        tempt=tempt+betas(j)*taf1(m,j)/(1+taf1(m,j))*pi(1,m,j);
    end
end
x(1)=x(1)/(1-tempt);

% tax_income = betas_ns(2:end,:).*(extax1./(1+extax1).*squeeze(pi(2:end,1,:))).*x_ns(2:end,:);
% x(1) = x(1)+sum(tax_income(:));

x_ns = x.*ones(1,ns);


for m=1:n % destination country
    for j=1:n
        pval0(j,:)= squeeze(xnimatrix(m,j,:)); % shipped from j to m
    end
    %pval1 = (sum(squeeze(xnimatrix(m,:,:)),2)).^(-betas./theta);
    pval1 = (sum(pval0,1)).^(-betas./theta);
    tp0 = cumprod(pval1);
    P(m,1) = tp0(end);
end

pi11 = (squeeze(pi(1,1,:)));

if (n==2)
    pim1= (squeeze(pi(2:n,1,:)))'; % squeeze(pi(1,1,:)), we got an column; squeeze pi(:,1,1) also column; squeeze(pi(1,:,1)) row
    exm = (squeeze(sextax(2:end,1,:)))';
    tafm = (squeeze(staf(1,2:end,:)))';
else
    pim1= (squeeze(pi(2:n,1,:)));
    exm = (squeeze(sextax(2:end,1,:)));
    tafm = (squeeze(staf(1,2:end,:)));
end

Lp(1,:) = theta/(1+theta)*betas./w(1).*(pi11'*x(1)+sum(1./(1+exm).*pim1.*x(2:end),1));

for i=2:n
    pi1n = (squeeze(pi(1,i,:)))';
    if (n==2)
        pimn= (squeeze(pi(2:end,i,:)))';
    else
        pimn= (squeeze(pi(2:end,i,:)));
    end
    Lp(i,:) = theta/(1+theta)*betas./w(i).*(pi1n*x(1)./(1+tafm(i-1,:))+sum(pimn.*x_ns(2:end,:),1));
end

Ls = Lp+rL;

uc = x(1)^(-sigma);

for i=1:n
    Mn(i,:) = w(i)./alpha(i,:);
end

%% construct matrix to solve 3*ns*(n-1)+n-2 with % gammaL1=0 and gamma4=uc
% gammaLn:ns*(n-1)-1
% gammarn:ns*(n-1)
% gammaTn: ns*(n-1)
% muN: n-1

% Note there are (n-1)*ns-1 gammaLnj, we will set the sector with
% CA has gammaLnj=0 since it for sure has Nj>0
% j is columun,
% inn: index on country
% in1: country index on matrix

tLs = Ls(2:end,:);
ttt=find(tLs==0);
gg = 1:(n-1)*ns;
gg(ttt)=[];
NG = length(gg);

% delete the multiplier gammaLnj with highest alpha for the second country
%if (n==2)
    [hs,iss]=max(alpha(2,:));
% else
% [hs,iss]=max(alpha(3,:));
% end
%ttt=find(gg==2+iss);
ttt=find(gg==iss);
cgg = gg;
cgg(ttt) = [];
cNG = length(cgg);


%gammaLn(cNG), gammarn(ns*(n-1)),gammaTn(ns*(n-1)), muN(n-1)
lhs_matrix = zeros(cNG+2*ns*(n-1)+n-1,cNG+2*ns*(n-1)+n-1);
rhs_matrix = zeros(cNG+2*ns*(n-1)+n-1,1);

%%  go through each row
%FOC on Tnj
for m=1:NG % FOC for Tnj
    j=gg(m);
    in1 = floor((j-1)/ns)+1;% country index on the matrix
    iss = j-(in1-1)*ns; % sector index
    inn = in1+1; % real country index
    
    %gammaLi
    pimn = squeeze(pi(:,inn,iss));
    for inp1 = 1:n-1 % country index on the matrix
        inp = inp1+1; % real country index
        
        mf = find(cgg==(inp1-1)*ns+iss);
        if (isempty(mf)==0)
            pimi = squeeze(pi(:,inp,iss));
            lhs_matrix(m,mf)=-betas(iss)*theta*sum(pimi(2:end).*pimn(2:end).*x(2:end));
        end
    end
    
    % gammaLn
    mf = find(cgg==(in1-1)*ns+iss);
    if (isempty(mf)==0)
        lhs_matrix(m,mf)=lhs_matrix(m,mf)+betas(iss)*theta*sum(pimn(2:end).*x(2:end));
    end
    
    % gammarn
    lhs_matrix(m,cNG+(in1-1)*ns+iss)=-w(inn)*Lp(inn,iss)/T(inn,iss);
    
    % gammari
    for inp1=1:n-1 % gammari index in the matrix
        for isp=1:ns
            inp = inp1+1;
            mf = cNG+(inp1-1)*ns+isp;
            lhs_matrix(m,mf)=lhs_matrix(m,mf)-beta*(1-delta)*(1-sigma)*Mn(inp,isp)*betas(iss)*pi(inp,inn,iss);
            % here we ignore dGi/dTni
        end
    end
    
    % gammaTn
    mf = cNG+(n-1)*ns+(in1-1)*ns+iss; % mun (n-1)
    lhs_matrix(m,mf)=-theta*T(inn,iss)+theta*beta*(1-delta)*T(inn,iss);

    % RHS
    pim1 = squeeze(pi(2:end,1,iss));
    pimn = squeeze(pi(2:end,inn,iss));
    taum = squeeze(sextax(2:end,1,iss));
    rhs_matrix(m,1) =-uc*betas(iss)*pi(1,inn,iss)*x(1)+uc*theta*betas(iss)*sum(taum./(1+taum).*pim1.*pimn.*x(2:end));
    %-rho*GTn(inn)/alpha(inn,iss)*(1+gL)*T(inn,iss)
end


%% FOC on Lnp
for m=1:(n-1)*ns
    
    in1 = floor((m-1)/ns)+1;% country index on the matrix
    iss = m-(in1-1)*ns; % sector
    inn = in1+1;% country index
    
    indd = NG+m;
    
    % gammaLn
    mf = find(cgg==m);
    if (isempty(mf)==0)
        lhs_matrix(indd,mf)= -1;
    end
    
    % gammaLi
    for inp1=1:n-1
        for isp=1:ns
            inp = inp1+1;
            mf=find(cgg==(inp1-1)*ns+isp);
            if (isempty(mf)==0)
                pini = squeeze(pi(inn,inp,isp));
                lhs_matrix(indd,mf) = lhs_matrix(indd,mf) + betas(isp)*pini;
            end
        end
    end
        
    % gammarn
    mf = cNG+(in1-1)*ns+iss;
    lhs_matrix(indd,mf) = 1/(1+theta)/T(inn,iss);
    
    %gammarn-k
    for isp=1:ns
        mf = cNG+(in1-1)*ns+isp;
        lhs_matrix(indd,mf) =lhs_matrix(indd,mf) +beta*(1-delta)*sigma*Mn(inn,isp)/x(inn);
    end
    
    % mun
    mf = cNG+2*(n-1)*ns+in1;
    lhs_matrix(indd,mf) = -theta/(1+theta)/w(inn);
    

    %RHS
    taum = squeeze(sextax(inn,1,:));
    pin1 = squeeze(pi(inn,1,:));
   
    rhs_matrix(indd,1) = -uc*sum(betas.*taum'./(1+taum').*pin1');
end

%% Lnr
for m=1:(n-1)*ns
    
    in1 = floor((m-1)/ns)+1;% country index on the matrix
    iss = m-(in1-1)*ns; % sector
    inn = in1+1;% country index
    
    indd = NG+(n-1)*ns+m;
    
    % gammaTn
    mf = cNG+(n-1)*ns+(in1-1)*ns+iss; % mun (n-1)
    lhs_matrix(indd,mf)=alpha(inn,iss);
    
    % mun
    mf = cNG+2*(n-1)*ns+in1;
    lhs_matrix(indd,mf) = -1;
    
    rhs_matrix(indd,1) = 0;
end

%% wn
if (n>2)
    for inn=3:n % country index
        in1 = inn-1; % country index in the matrix 
        
        indd = NG+2*(n-1)*ns+(inn-2);
        
        % gammaLn
        for isp=1:ns
            ind2 =(in1-1)*ns+isp; %CHECK!!!
            mf=find(ind2==cgg);
            if (isempty(mf)==0)
                if (n==2)
                    pimn = squeeze(pi(2:end,inn,isp))';
                    
                else
                    pimn = squeeze(pi(2:end,inn,isp));
                    
                end
                lhs_matrix(indd,mf) = -(1+theta)/theta*w(inn)*Lp(inn,isp)-betas(isp)*theta*sum(pimn.*x(2:end));
            end
        end
        
        
        % gammaLi
        for inp1=1:n-1
            for isp=1:ns
                inp = inp1+1;
                mf=find(cgg==(inp1-1)*ns+isp);
                if (isempty(mf)==0)
                    
                    if (n==2)
                        pim1 = squeeze(pi(2:end,1,isp))';
                        pimi = squeeze(pi(2:end,inp,isp))';
                        pimn = squeeze(pi(2:end,inn,isp))';
                        
                    else
                        pim1 = squeeze(pi(2:end,1,isp));
                        pimi = squeeze(pi(2:end,inp,isp));
                        pimn = squeeze(pi(2:end,inn,isp));
                        
                    end
                    pini = squeeze(pi(inn,inp,isp));
                    
                    lhs_matrix(indd,mf) =lhs_matrix(indd,mf)+betas(isp)*theta*sum(pimi.*pimn.*x(2:end))...
                        +betas(isp)*pini*x(inn);
                end
            end
        end
       
        % gammarn
        for isp=1:ns
            mf = cNG+(in1-1)*ns+isp;
            lhs_matrix(indd,mf) = beta*(1-delta)*(sigma-1)*Mn(inn,isp);
        end
        
        %gammari
        for inp1=1:n-1
            for isp=1:ns
                inp = inp1+1;
                mf = cNG+(inp1-1)*ns+isp;
                piin = squeeze(pi(inp,inn,:));
                lhs_matrix(indd,mf)=lhs_matrix(indd,mf)+beta*(1-delta)*(1-sigma)*Mn(inp,isp)*sum(betas.*piin');
            end
        end
        
        pi1n = squeeze(pi(1,inn,:));
        if (n==2)
            pim1 = squeeze(pi(2:end,1,:))';
            pimn = squeeze(pi(2:end,inn,:))';
            taum = squeeze(sextax(2:end,1,:))';
        else
            pim1 = squeeze(pi(2:end,1,:));
            pimn = squeeze(pi(2:end,inn,:));
            taum = squeeze(sextax(2:end,1,:));
        end
        
        pin1 = squeeze(pi(inn,1,:));
        taun = squeeze(sextax(inn,1,:));
        % RHS

        rhs_matrix(indd,1)= uc*x(1)*sum(betas'.*pi1n)-uc*(sum(sum(betas_ns(2:end,:).*taum./(1+taum)*theta.*pim1.*pimn.*x_ns(2:end,:)))...
            +sum(betas'.*taun./(1+taun).*pin1.*x(inn)));
        
    end
end

%%
tgamma = lhs_matrix\rhs_matrix;
tgammaLn = zeros((n-1)*ns,1);
tgammaLn(cgg) = tgamma(1:cNG);
tgammarn = tgamma(cNG+1:cNG+(n-1)*ns);
tgammaTn = tgamma(cNG+(n-1)*ns+1:cNG+2*(n-1)*ns);
mun1 = tgamma(cNG+2*(n-1)*ns+1:cNG+2*(n-1)*ns+n-1);
gammaL1 = 0;
gamma4 = uc;

gammaL(1,1:ns) = gammaL1;
gammaL(2:n,1:ns) = reshape(tgammaLn,ns,n-1)';
gammarn(1,1:ns) = 0;
gammarn(2:n,1:ns) = reshape(tgammarn,ns,n-1)';
mun(1) = (gamma4-gammaL1)*(1+theta)/theta*w(1);
mun(2:n) = mun1;

gammaTn(1,1:ns) = 0;
gammaTn(2:n,1:ns) = reshape(tgammaTn,ns,n-1)';


%%

fval(1,1) = P(1)-1;
for i=2:n
    fval(i,1) = sum(rL(i,:)+Lp(i,:))-L(i);
end

fval(n+1:n+ns,1) =rL(1,:)./Ls(1,:)-rbest;

for i=1:n-1
    inn = i+1;
%     fval(n+ns+(i-1)*ns+1:n+ns+i*ns,1)= w(inn)./alpha(inn,:)'- (1/theta)*w(inn).*(Lp(inn,:)')./(T(inn,:)'*(1+gL))...
%         -rho*w(inn)./alpha(inn,:)';
    
    fval(n+ns+(i-1)*ns+1:n+ns+i*ns,1) =rL(inn,:)./Ls(inn,:)-rbest;
end

pi11 = squeeze(pi(1,1,:));
pi1n = squeeze(pi(1,:,:));

for i=1:n-1
    inn = i+1;
    
    taum = squeeze(sextax(inn,1,:));
    pin1 = squeeze(pi(inn,1,:));
    pini = squeeze(pi(inn,2:end,:));
    ins = n+ns+(n-1)*ns+(i-1)*ns;
    
    tafm = squeeze(staf(1,inn,:));
    
    fval(ins+1:ins+ns,1)=(1+tafm)-(gamma4-gammaL(inn,:)')/gamma4;% tariff
    
    
    tempt=theta*sum(gammaL(2:end,:).*pini,1)';
    tempt2 = beta*(1-delta)*(sigma-1)*sum(gammarn(inn,:).*Mn(inn,:))/x(inn);
    
    fval(ns+ins+1:ns+ins+ns,1)=(1+taum)-uc*(1+theta*(1-pin1))./(uc*theta*(1-pin1)-tempt+tempt2);
end







